---
title: "10_combinedresults"
author: "Ye Shen, Radhika Tampi, Raj Vatsa"
date: "12/7/2021"
output: pdf_document
---
This rmarkdown file shows code to prepare combined figures and tables from the extension analyses for use in the final version of the paper. 

```{r}
library(tidyverse)
library(patchwork)
library(ggplot2)
```

# Combined Figure for Children
## Read in ggplot objects from each extension analysis for children
```{r}
# Matching panel
fig_child_match <-readRDS("table_fig/child/fig_child_match.rds") +
                  theme(legend.position ="None", 
                        plot.title=element_blank(),
                        axis.title.x=element_blank())+
                  labs(subtitle="Matching on Household Wealth Index")

# Weighting panel
fig_child_wt<-readRDS("table_fig/child/fig_child_wt.rds")+
                  theme(legend.position ="None", 
                        plot.title=element_blank(),
                        axis.title.x=element_blank())+
                  labs(subtitle="Balance Weighting on Household Wealth Index")

# Mixed effects panel
fig_child_mixed<-readRDS("table_fig/child/fig_child_mixed.rds")+
                  theme(plot.title=element_blank())+
                  labs(subtitle="Mixed Effects DiD with Village-Level Random Effects")

# Original panel
final_table2_children_mixed <- 
  readRDS("table_fig/child/final_table2children_mixed.rds")

num_var <- 1
fig_child_orig <- as.data.frame(matrix(nrow = num_var * 4, ncol = 4))
fig_child_orig$V1 <- round(c(
  as.numeric(gsub(" .*$", "", final_table2_children_mixed$m_2014)), 
  as.numeric(gsub(" .*$", "", final_table2_children_mixed$m_2016))), 2)
fig_child_orig$V2 <- rep(c("Intervention", "Non-Intervention"), num_var * 2)
fig_child_orig$V3 <- rep(c("2014 survey", "2016 survey"), each = num_var * 2)
fig_child_orig$V4 <- 
  rep("Any Sick-Child Visits to a PHF or CHW", num_var * 4) ### RENAME BASED ON OUTCOME ###
colnames(fig_child_orig) <- c("mean", "group", "year", "indicator")
fig_child_orig$group <- factor(fig_child_orig$group, levels = c("Non-Intervention", "Intervention"))

fig_child_orig <- 
  ggplot(fig_child_orig, 
         aes(x = mean, y = group, group = interaction(group, indicator))) +
  geom_point(aes(colour = year), size = 3.5) +
  scale_colour_manual(breaks = fig_child_orig$year, 
                      values = c("2014 survey" = "black", "2016 survey" = "red")) +
  geom_line() +
  labs(
    x = "", y = "",
    subtitle="Ezran et al. (2019): Unadjusted and Adjusted for Household Wealth"
  ) +
  theme_bw() +
  scale_x_continuous(limits = c(0, 100), breaks = seq(0, 100, 10), expand = c(0, 0)) +
  scale_y_discrete(position = "right") +
  theme(
    legend.position = "None", axis.text = element_text(size = 10), 
    panel.grid.minor = element_blank(), legend.title = element_blank(),
    title = element_text(size = 12, face = "bold"), 
    axis.title = element_text(size = 13),
    legend.text = element_text(size = 12)
  )+
  annotate("text", label="}", x=70, y=1.5, size =20, fontface="plain")+
  annotate("text", label=paste0("Unadjusted DiD (p-val): \n",
                                final_table2_children_mixed$orig_DiD[2]), x=85, y=1.65) + 
  annotate("text", label=paste0("Adjusted DiD (p-val): \n",
                                final_table2_children_mixed$orig_DiD_adj[2]), x=85, y=1.15)
```
## Combine children figures
```{r}
fig_combined_child<-fig_child_orig/fig_child_match/fig_child_wt/fig_child_mixed + 
  plot_annotation(
    title="Figure 4. Any Sick-Child Visits to a PHF or CHW",
    tag_levels = 'A',
    caption="Note: Adjusted DiD included time-varying household wealth index as a control"
        ) &
    theme(plot.tag = element_text(size = 10), text = element_text( face="bold"))
ggsave("table_fig/paper/fig4_combined_child.png",plot=fig_combined_child,width = 200,height=300, units="mm")
  
```

# Combined Figure for Antenatal Services
## Read in ggplot objects from each extension analysis for pregnant women
```{r}
# Matching panel
fig_ante_match<-readRDS("table_fig/antenatal/fig_ante_match.rds") +
                  theme(legend.position ="None",
                        plot.title=element_blank(),
                        axis.title.x=element_blank())+
                  labs(subtitle="Matching on Household Wealth Index")

# Weighting panel
fig_ante_wt<-readRDS("table_fig/antenatal/fig_ante_wt.rds")+
                  theme(legend.position ="None", 
                        plot.title=element_blank(),
                        axis.title.x=element_blank())+
                  labs(subtitle="Balance Weighting on Household Wealth Index")

# Mixed effects panel
fig_ante_mixed<-readRDS("table_fig/antenatal/fig_antenatal_mixed.rds")+
                  theme(plot.title=element_blank())+
                  labs(subtitle="Mixed Effects DiD with Village-Level Random Effects")

# Original panel
final_table3_antenatal_mixed <- 
  readRDS("table_fig/antenatal/final_table3antenatal_mixed.rds")

num_var <- 1
fig_ante_orig <- as.data.frame(matrix(nrow = num_var * 4, ncol = 4))
fig_ante_orig$V1 <- round(c(
  as.numeric(gsub(" .*$", "", final_table3_antenatal_mixed$m_2014)), 
  as.numeric(gsub(" .*$", "", final_table3_antenatal_mixed$m_2016))), 2)
fig_ante_orig$V2 <- rep(c("Intervention", "Non-Intervention"), num_var * 2)
fig_ante_orig$V3 <- rep(c("2014 survey", "2016 survey"), each = num_var * 2)
fig_ante_orig$V4 <- 
  rep("Any Antenatal Care Visits to a PHF or CHW", num_var * 4) ### RENAME BASED ON OUTCOME ###
colnames(fig_ante_orig) <- c("mean", "group", "year", "indicator")
fig_ante_orig$group <- factor(fig_ante_orig$group, levels = c("Non-Intervention", "Intervention"))

fig_ante_orig <- 
  ggplot(fig_ante_orig, 
         aes(x = mean, y = group, group = interaction(group, indicator))) +
  geom_point(aes(colour = year), size = 3.5) +
  scale_colour_manual(breaks = fig_child_orig$year, 
                      values = c("2014 survey" = "black", "2016 survey" = "red")) +
  geom_line() +
  labs(
    x = "", y = "",
    subtitle="Ezran et al. (2019): Unadjusted and Adjusted for Household Wealth"
  ) +
  theme_bw() +
  scale_x_continuous(limits = c(0, 100), breaks = seq(0, 100, 10), expand = c(0, 0)) +
  scale_y_discrete(position = "right") +
  theme(
    legend.position = "None", axis.text = element_text(size = 10), 
    panel.grid.minor = element_blank(), legend.title = element_blank(),
    title = element_text(size = 12, face = "bold"), 
    axis.title = element_text(size = 13),
    legend.text = element_text(size = 12)
  )+
  annotate("text", label="{", x=50, y=1.5, size =20, fontface="plain")+
  annotate("text", label=paste0("Unadjusted DiD (p-val): \n",
                                final_table3_antenatal_mixed$orig_DiD[2]), x=35, y=1.65) + 
  annotate("text", label=paste0("Adjusted DiD (p-val): \n",
                                final_table3_antenatal_mixed$orig_DiD_adj[2]), x=35, y=1.15)
```
## Combine antenatal figures
```{r}
fig_combined_ante<-fig_ante_orig/fig_ante_match/fig_ante_wt/fig_ante_mixed + 
  plot_annotation(
    title="Figure 5. Any Antenatal Care Visits to a PHF or CHW",
    tag_levels = 'A',
    caption="Note: Adjusted DiD included time-varying household wealth index as a control"
        ) &
    theme(plot.tag = element_text(size = 10), text = element_text( face="bold"))
ggsave("table_fig/paper/fig5_combined_ante.png",plot=fig_combined_ante,width = 200,height=300, units="mm")
```

# Combined Figure for Perinatal Services
## Read in ggplot objects from each extension analysis for perinatal services
```{r}
# Matching panel
fig_peri_match<-readRDS("table_fig/perinatal/fig_peri_match.rds") +
                  theme(legend.position ="None",
                        plot.title=element_blank(),
                        axis.title.x=element_blank())+
                  labs(subtitle="Matching on Household Wealth Index")

# Weighting panel
fig_peri_wt<-readRDS("table_fig/perinatal/fig_peri_wt.rds")+
                  theme(legend.position ="None", 
                        plot.title=element_blank(),
                        axis.title.x=element_blank())+
                  labs(subtitle="Balance Weighting on Household Wealth Index")

# Mixed effects panel
fig_peri_mixed<-readRDS("table_fig/perinatal/fig_perinatal_mixed.rds")+
                  theme(plot.title=element_blank())+
                  labs(subtitle="Mixed Effects DiD with Village-Level Random Effects")

# Original panel
final_table3_perinatal_mixed <- 
  readRDS("table_fig/perinatal/final_table3perinatal_mixed.rds")

num_var <- 1
fig_peri_orig <- as.data.frame(matrix(nrow = num_var * 4, ncol = 4))
fig_peri_orig$V1 <- round(c(
  as.numeric(gsub(" .*$", "", final_table3_perinatal_mixed$m_2014)), 
  as.numeric(gsub(" .*$", "", final_table3_perinatal_mixed$m_2016))), 2)
fig_peri_orig$V2 <- rep(c("Intervention", "Non-Intervention"), num_var * 2)
fig_peri_orig$V3 <- rep(c("2014 survey", "2016 survey"), each = num_var * 2)
fig_peri_orig$V4 <- 
  rep("Deliver at a PHF", num_var * 4) ### RENAME BASED ON OUTCOME ###
colnames(fig_peri_orig) <- c("mean", "group", "year", "indicator")
fig_peri_orig$group <- factor(fig_peri_orig$group, levels = c("Non-Intervention", "Intervention"))

fig_peri_orig <- 
  ggplot(fig_peri_orig, 
         aes(x = mean, y = group, group = interaction(group, indicator))) +
  geom_point(aes(colour = year), size = 3.5) +
  scale_colour_manual(breaks = fig_child_orig$year, 
                      values = c("2014 survey" = "black", "2016 survey" = "red")) +
  geom_line() +
  labs(
    x = "", y = "",
    subtitle="Ezran et al. (2019): Unadjusted and Adjusted for Household Wealth"
  ) +
  theme_bw() +
  scale_x_continuous(limits = c(0, 100), breaks = seq(0, 100, 10), expand = c(0, 0)) +
  scale_y_discrete(position = "right") +
  theme(
    legend.position = "None", axis.text = element_text(size = 10), 
    panel.grid.minor = element_blank(), legend.title = element_blank(),
    title = element_text(size = 12, face = "bold"), 
    axis.title = element_text(size = 13),
    legend.text = element_text(size = 12)
  )+
  annotate("text", label="}", x=45, y=1.5, size =20, fontface="plain")+
  annotate("text", label=paste0("Unadjusted DiD (p-val): \n",
                                final_table3_perinatal_mixed$orig_DiD[2]), x=65, y=1.65) + 
  annotate("text", label=paste0("Adjusted DiD (p-val): \n",
                                final_table3_perinatal_mixed$orig_DiD_adj[2]), x=65, y=1.15)
```

## Combine perinatal figures
```{r}
fig_combined_peri<-fig_peri_orig/fig_peri_match/fig_peri_wt/fig_peri_mixed + 
  plot_annotation(
    title="Figure 6. Delivery at a PHF",
    tag_levels = 'A',
    caption="Note: Adjusted DiD included time-varying household wealth index as a control"
        ) &
    theme(plot.tag = element_text(size = 10), text = element_text( face="bold"))
ggsave("table_fig/paper/fig6_combined_peri.png",plot=fig_combined_peri,width = 200,height=300, units="mm")


```

## Combine village-level variation figures

```{r}
fig_min_village_child <- readRDS("table_fig/child/village_min_child_plot.rds") + 
  theme(legend.position ="None",
        axis.title.x=element_blank())
fig_median_village_child <- readRDS("table_fig/child/village_median_child_plot.rds") + 
  theme(legend.position ="None",
        plot.title=element_blank(),
        axis.title.x=element_blank())
fig_max_village_child <- readRDS("table_fig/child/village_max_child_plot.rds") + 
  theme(legend.position ="None",
        plot.title=element_blank(),
        axis.title.x=element_blank())

fig_min_village_antenatal <- readRDS("table_fig/antenatal/village_min_antenatal_plot.rds") + 
  theme(legend.position ="None",
        axis.title.x=element_blank())
fig_median_village_antenatal <- readRDS("table_fig/antenatal/village_median_antenatal_plot.rds") + 
  theme(legend.position ="None",
        plot.title=element_blank(),
        axis.title.x=element_blank())
fig_max_village_antenatal <- readRDS("table_fig/antenatal/village_max_antenatal_plot.rds") + 
  theme(legend.position ="None",
        plot.title=element_blank(),
        axis.title.x=element_blank())

fig_min_village_perinatal <- readRDS("table_fig/perinatal/village_min_perinatal_plot.rds")
fig_median_village_perinatal <- readRDS("table_fig/perinatal/village_median_perinatal_plot.rds") + 
  theme(plot.title=element_blank(),
        legend.position ="None")
fig_max_village_perinatal <- readRDS("table_fig/perinatal/village_max_perinatal_plot.rds") + 
  theme(plot.title=element_blank(),
        legend.position ="None")

fig_combined_variation <- (fig_min_village_child + fig_median_village_child + fig_max_village_child) / 
  (fig_min_village_antenatal + fig_median_village_antenatal + fig_max_village_antenatal) / 
  (fig_min_village_perinatal + fig_median_village_perinatal + fig_max_village_perinatal) +  
  plot_annotation(
    title="Figure 7. Variation in Estimated Intervention Effect Across Villages",
    tag_levels = 'A',
    caption="Note: DiD estimate includes time-varying household wealth index as a control"
        ) &
    theme(plot.tag = element_text(size = 10), text = element_text( face="bold"))
ggsave("table_fig/paper/fig7_combined_mixed_effect_variation.png",
       plot = fig_combined_variation, width = 300, height=300, units="mm")
```

